# 23_utas2003.R
# Purpose: Ideological Extremism and Political Participation in Japan
# Created: 2021-3-14 Taka-aki Asano
# Last Modified: 2021-10-10

# package
require("ltm")
require("plink")
require("haven")
require("dplyr")
require("tidyr")
require("psych")


# Voters
## read dataset
voter2003 <- read_sav(
  "utas060104.sav", 
  encoding = "CP932", 
  user_na = FALSE
)

## specify variables
variables2003 <- c("q021501", "q021502", "q021503", "q021504", "q021505",
                   "q021506", "q021507", "q021508", "q021509", "q021510", 
                   "q021511", "q021512", "q021513")

## handle missing data
for(i in variables2003) {
  n <- voter2003[,i] == 1.5 | voter2003[,i] == 2.5 | 
    voter2003[,i] == 3.5 | voter2003[,i] == 4.5
  n[is.na(n)] <- FALSE
  if (sum(n) != 0) {
    voter2003[n, i] <- NA
  }
}
voter2003 <- voter2003[rowSums(voter2003[,variables2003], na.rm = TRUE) != 0,]

## estimate item parameters
irt2003_voter <- ltm::grm(voter2003[,variables2003], IRT.param = FALSE)
irt2003_voter ## viewing

## estimate voters ideology
score2003 <- ltm::factor.scores(irt2003_voter, voter2003[,variables2003])


# rescale
## common item
common2009 <- c("Q013802", "Q013804", "Q013805", "Q013809", 
                "Q013810", "Q013811", "Q013812", "Q013816")
common2003 <- c("q021501", "q021503", "q021504", "q021506", 
                "q021507", "q021508", "q021509", "q021513")
index <- data.frame(
  group1 = match(common2009, colnames(voter2009[,variables2009])), 
  group2 = match(common2003, colnames(voter2003[,variables2003]))
)
index

## list of parameters
### 2009
res2009 <- coef(irt2009_voter)
res2009 <- res2009[,c(5,1:4)]
colnames(res2009) <- c("a","b1","b2","b3","b4")
### 2003
res2003 <- coef(irt2003_voter)
res2003 <- res2003[,c(5,1:4)]
colnames(res2003) <- c("a","b1","b2","b3","b4")
### merge
pm <- list(res2009, res2003)

## list of scale
rescat <- list(rep(5, nrow(res2009)), 
               rep(5, nrow(res2003)))

## as.irt.pars
pm2009 <- as.poly.mod(n = nrow(res2009), model = "grm", 
                      items = 1:nrow(res2009))
pm2003 <- as.poly.mod(n = nrow(res2003), model = "grm", 
                      items = 1:nrow(res2003))
res <- as.irt.pars(pm, common = index, cat = rescat, 
                   poly.mod = list(pm2009, pm2003), 
                   location = FALSE)

## list of voters' position
ideology <- list(score2009$score.dat$z1, score2003$score.dat$z1)

## rescale
plink_out <- plink(res, common = index, base.grp = 1, 
                   rescale = "MS", ability = ideology)
summary(plink_out)

## estimate voters ideology (rescale)
score2003_plink <- link.ability(plink_out)$group2
voter2003$Ideology <- -1 * score2003_plink

## Median
median(voter2003$Ideology)
by(voter2003$Ideology, voter2003$q020600, describe)
